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ABSTRACT 

We present general relativistic magnetohydrodynamic (GRMHD) numerical simu- 
lations of the accretion flow around the supermassive black hole in the Galactic centre, 
Sagittarius A* (Sgr A*). The simulations include for the first time radiative cooling 
processes (synchrotron, bremsstrahlung, and inverse Compton) self-consistently in the 
dynamics, allowing us to test the common simplification of ignoring all cooling losses 
in the modeling of Sgr A*. We confirm that for Sgr A*, neglecting the cooling losses is 
a reasonable approximation if the Galactic centre is accreting below ~ 1O -8 M0 yr _1 
i.e. M < 10~ 7 M-Edd- But above this limit, we show that radiative losses should be 
taken into account as significant differences appear in the dynamics and the result- 
ing spectra when comparing simulations with and without cooling. This limit implies 
that most nearby low-luminosity active galactic nuclei are in the regime where cooling 
should be taken into account. 

We further make a parameter study of axisymmetric gas accretion around the 
supermassive black hole at the Galactic centre. This approach allows us to investigate 
the physics of gas accretion in general, while confronting our results with the well 
studied and observed source, Sgr A*, as a test case. We confirm that the nature of 
the accretion flow and outflow is strongly dependent on the initial geometry of the 
magnetic field. For example, we find it difficult, even with very high spins, to generate 
powerful outflows from discs threaded with multiple, separate poloidal field loops. 

Key words: accretion discs - black hole physics - MHD - radiation mechanisms: 
thermal - relativistic processes - methods: numerical - Galaxy: center - Galaxy: nu- 
cleus - galaxies: jets. 



1 INTRODUCTION 

Super-massive black holes (SMBHs) of millions to billions 
of solar masses are believed to exist in the centre of most 
galaxies. The Galactic center black hole candidate, Sgr A*, 
is the closest and best studied SMBH, making it the perfect 
source to test our understanding of galactic nuclei systems 
in general. This compact object was first observed as a radio 



cretion rate has also been constrained by polarisation mea- 



source by Balick & Brown (19741. Since then, observations 
have constrained important parameters of Sgr A*, such as 
its mass and distance, estimated at M = 4.3 ± 0.5 x 1O 6 M0 
and D = 8.3 ± 0.35 kpc, respectively ( Reidl 1993 Schodel 



et al.|2002| |Ghez et al.|[2008] |Gillessen et al.|2009| >. The ac 
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surements, using Faraday rotation arguments ( Aitken et al 



2000 Bower et al. 2003 



Marrone et al. 20071 and is esti- 
a < M < 2 x 1O" 7 M yr _1 . 

inclination, and 



mated to be in the range 2 x 10 
Other key parameters, such as the spin 
magnetic field configuration, are still under investigation. 

Multi-wavelength observations of Sgr A* have been 
performed from the radio to the gamma ray (see reviews 
by |Melia fc Falckel [2001] |Genzel et al~| |20l0l and refer- 
ences therein). More recently, important progress has been 
achieved in the infrared (IR; e.g., Schoedel et al.pOll l and 
sub-millimeter (sub-mm) domains. All observations agree 
that Sgr A* is a very under-luminous and weakly accret- 
ing black hole; indeed its accretion rate is lower than has 
been observed in any other accreting system. 

In the near future, the next milestone will be observ- 
ing the first black hole shadow from Sgr A* (Falcke et al. 
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2000 



Dexter et al.|2010 1 with the proposed "Event Horizon 



Telescope" dDoeleman et al. 1120081 120091 iFish et al.||2011l 



thanks to the capabilities of very long base interferometry 
at sub-mm wavelengths. Such a detection would be the first 
direct evidence for a black hole event horizon, and may also 
constrain the spin of Sgr A*. In order to make accurate 
predictions for testing with the Event Horizon Telescope, 
however, we need to have reliable models for the plasma 
conditions and geometry in the accretion (in/out)flow. Gen- 
eral relativistic magnetohydrodynamic (MHD) simulations 
offer significant promise for this class of study, as they can 
provide both geometrical and spectral predictions. 

Sgr A* has already been modeled in several numerical 



studies (e.g 


.,|Goldston et al. 2005 


Moscibrodzka et al.|2009 


Dexter et al.|2009 2010 Hilburn et al. 2010 Shcherbakov| 


et al.j 2010 


Shiokawa et al.|2012 


?■ 


Dexter & Fragile 2012 


1. 



All of these models consist of two separate codes: a GRMHD 
code describing the dynamics, and then a subsequent code 
to calculate the radiative emission based on the output of 
the first one. The under-luminous and under-accreting state 
of Sgr A* (Lboi ^ 10" 9 L E dd and L X -R ay ^ 10 33 erg/s in 
the 0.5 to 10 keV band, where -Z^Edd is the Eddington lu- 
minosity; Baganoff et al. 20031, is the common argument 



given to justify ignoring the cooling losses and simplify the 
description of Sgr A*. Even though this approach seems rea- 
sonable, especially for the peculiar case of Sgr A*, we pro- 
pose to quantify to what extent it really applies. This ques- 
tion is especially important if one wants to extend these 
studies to more typical nearby low-luminosity active galac- 
tic nuclei (LLAGN). An important example is the nuclear 
black hole in M87, which is the other major target for 
mm-VLBI and is significantly more luminous than Sgr A* 



10" 



10" 4 ). 



To test whether or not ignoring the radiative cooling 
losses is a reasonable approximation, we need to take into 
account the radiative losses self-consistently in the dynam- 
ics, which is now possible with the Cosmos++ astrophysical 



fluid dynamics code (Anninos et al. 2005 Fragile & Meier 



20091. By allowing the gas to cool, energy can be liberated 
from the accretion flow, potentially changing the dynamics 
of the system. The basic result that we show in this paper 
is that cooling is indeed negligible at the lowest end of the 
possible accretion rate range of Sgr A*, but plays an in- 
creasingly important role in the dynamics and the resulting 
spectra for higher accretion rates (relevant for most nearby 
LLAGN and even at the high end of the range of Sgr A*). 

This paper is organized as follows: In Section 2 we de- 
scribe the new version of Cosmos++ used to perform this 
study. In Section 3 we present the initial set-up of the simu- 
lations. In Section 4 we show the importance of the accretion 
rate parameter and assess the effect of including radiative 
cooling. In Section 5 we perform a parameter survey to in- 
vestigate the influence of the initial magnetic field configu- 
ration and the spin of the SMBH (we also discuss the effect 
of a retrograde spin on the dynamics) on the resulting ac- 
cretion disc structure and outflow. In Section 6 we end with 
our conclusions and outlook. The spectra generated from 
these simulations are analyzed in detail in a companion pa- 
per (Drappeau et al., in preparation, hereafter referred to in 
the text as D12). 



2 NUMERICAL METHOD 
2.1 GRMHD Equations 

Within Cosmos++ we solve the following set of conserva- 
tive, general relativistic MHD equations, including radiative 
cooling, 



OtD + d^DV 1 ) 

dtSj + di(yf^g Tj) 
d t B j + di{B j V z -B l V 3 ) 



= 







_ _ (1) 
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where D = Wp is the generalized fluid density, W = y/—gu° 
is the generalized boost factor, V 1 = u' /u° is the transport 
velocity, u M = g^u^ is the fluid 4-velocity, is the 4- 
metric, g is the 4-metric determinant, 

s = -V^gTg 

= -Wu (ph + 2P B )- (P + P B ) + ^gb°b (5) 

is the total energy density, h — 1 + e + Pj p is the specific 
enthalpy, e is the specific internal energy, P is the fluid pres- 
sure, Pb is the magnetic pressure, T£ is the stress-energy 
tensor, A is the cooling function, and 

r,0 



-cfT} = W Uj {ph + 2P B )~ V^gFbj (6) 

is the covariant momentum density. With indices, V indi- 
cates the geometric connection coefficients of the metric. 
Without indices, F is the adiabatic index. For this work we 
use an ideal gas equation of state (EOS), 



P = (T - l)pe , 



(7) 



with F = 5/3. 

There are multiple representations of the magnetic field 
in our equations: b M is the magnetic field measured by an 
observer comoving with the fluid, which can be defined in 
terms of the dual of the Faraday tensor = u v * F^ v , and 
B 3 = \J—gP> 3 is the boosted magnetic field 3-vector. The 
un-boosted magnetic field 3-vector B % = * is related to 
the comoving field by 

B l = u°b* - u l b° . (8) 

The magnetic pressure is Pb = b 2 /2 = 6 M 6 M /2. Note that, 
unlike in previous versions of Cosmos++, we have absorbed 
the factor of y4n into the definition of the magnetic fields, 
so-called Lorentz-Heaviside units (£?lh = -Bcgs/V^Lr). 



2.2 GRMHD Solver 

We note that the MHD equations as written are all in the 
form of conservation equations 



au(p) + ftF'(p) = s(p) 



(9) 



where U, F 1 , and S represent the conserved quantities, 
fluxes, and source terms, respectively. These are solved using 
a new High Resolution Shock Capturing (HRSC) scheme re- 
cently added to the Cosmos++ computational astrophysics 



code (Anninos et al. 20051. The new HRSC scheme is 



modeled after the original non-oscillatory central difference 



(NOCD) scheme of Cosmos ( |Anninos fcFr agilc 200j|. It 
also has many of the same elements as the publicly available 
HARM code ( jGammie et~aT1[2003} [Noble et al.||2006[ ), al- 
though with staggered magnetic fields (instead of HARM's 
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zone-centered fields) and the inclusion of a self-consistent, 
physical cooling model. More details of the new scheme are 
provided in a separate paper ( Fragile et al.|2012 |. 

Briefly, as in other, similar conservative codes, the 
fluxes, F l , are determined using the Harten-Lax-Van Leer 



modifications to the original Esin et al. ( 1996 1 cooling func- 



tion are described in Fragile & Meier ( 2009 1 and are also 



(HLL) approximate Riemann solver (Harten et al.|[l983 1 



p CminF R -\- C max F L CmaxCminftJ R U l) (10) 

' rn a x 

A slope-limited parabolic extrapolation gives Pr and P_l, 
the primitive variables at the right- and left-hand side of 
each zone interface. From P_r and Pl , we calculate the right- 
and left-hand conserved quantities (\Jr and U_l), the fluxes 
Fr — F(Pr) and Fl = F(Pj,), and the maximum right- 
and left-going waves speeds, c±,r and c±,l. The bounding 
wave speeds are then c max = max(0, c +> r, c +i l) and c m i n = 
-min(0, c_,ij, c_,l). 

For stability, the equations are integrated using a stag- 
gered leapfrog method (2nd order in time). First a half-time 
step, At/2, is taken to project the conserved variables U n 
forward to n + 1/2. From these, a new set of primitives 
pn+1/2 can k e corn p U ted. These intermediate primitives are 
then used in calculating the fluxes and source terms needed 
for the full time step, At. To find the new primitive vari- 
ables, P, after each update cycle, we use either the 2D or 
lDw numerical methods of Noble et al. (20061. 



Along with the evolution equations, we must also sat- 
isfy the divergence- free constraint djB J — 0. To accomplish 
this, Cosmos++ now has the option to use a staggered 
magnetic field with a Constrained-Transport (CT) update 
scheme akin to the Newtonian version described in lStone~fel 



Gardiner (2009) 



2.3 Cooling 

In the present work, we assume the whole system is op- 
tically thin, i.e. radiation escapes freely from the system. 
However, for the calculation of the radiation we consider the 
appropriate optical depth of the gas at a given location and 
time, which depends on the state. This approximation takes 
into account the optical depth without performing radiative 
transfer, and is valid as long as the (assumed thermal) peak 
of the radiating particle distribution corresponds to energies 
greater than the self-absorption frequency, which is almost 
always the case for the regions under study. 

The radiative cooling term, A, in equations (2| and Q, 
is the cooling function introduced to Cosmos++ in fragile 
fc Meier| ( |2009[ ) and is based on |Esin et al] ( |1996| >, which in- 
cludes treatments of bremsstrahlung, synchrotron, and the 
inverse-Compton enhancement of each of these. Since Sgr 
A* is optically thin, we do not need to worry about multi- 
scattering events in the treatment of the inverse Compton 
emission. The cooling function is calculated locally in each 
zone from the fluid density, electron temperature, magnetic 
field strength, and scale height [A = A(p, T e , b 2 , H)]. The lo- 
cal electron temperature is recovered by assuming a constant 
ion-to-electron temperature ratio, T;/T e , (fixed for each sim- 
ulation) and calculating the ion temperature from the fluid 
variables Ti = pmnP/kBp, where p, = 1.69 is the mean 
molecular weight (appropriate for Solar abundances), mn is 
the mass of hydrogen, and ks is Boltzman's constant. The 
temperature scale height is defined as H = T* /|V(T*)|. Our 



described in more detail in our companion paper (D12). 

In taking the ion and electron temperatures to be held 
in a fixed ratio, we are assuming there is some efficient 
coupling process at work between the ions and electrons. 



The same was done in Esin et al. ( 1996 1 and Fragile & 



Meier ( 2009 I , except here we relax the assumption some- 



what so that the coupling does not have to be exact (Ti/T e 
can be greater than 1). This procedure is not entirely self- 
consistent since the expression for ion-electron collisions in 



bremsstrahlung cooling [Equation (17) of Fragile & Meier 
(2009])] assumes the ions and electrons have the same tem- 
perature. Therefore, whenever T > T e , we are actually un- 
derestimating the amount of bremsstrahlung cooling. How- 
ever, because bremsstrahlung is such a minor contributor to 
the cooling, this omission is not significant in the context 
of our current work. Obviously, assuming a fixed ratio of 
Ti/T e is not likely to be physically accurate. However, it is 



the current standard in most simulations (e.g. Moscibrodzka 
et al. |2009" |Dexter et al.|20i0 I and recent comparisons with 
more sophisticated treatments has so far not shown large 
discrepancies (Dexter, Quataert, priv. comm.). 

Because the cooling time of the gas, i coo i = pe/A, can 
sometimes be shorter than the MHD timestep, A^mhd ~ 
Ax/V, where Ax and V are the characteristic zone length 
and velocity, respectively, we allow the cooling routine to op- 
erate on its own shorter timestep (subcycle), if necessary. To 
prevent runaway loops inside the cooling function, we limit 
the cooling subcycle to 4 steps, with a maximum change 
to the specific internal energy each subcycle of Ae = 0.5e. 
Even with these restrictions, the cooling has the potential 
to decrease the internal energy (and hence temperature) by 
up to nearly 95% each MHD cycle. 



2.4 Spectra 

Our method for generating spectra is described in detail 
in our companion paper (D12). However, we want to men- 
tion one important point related to how we present spectra 
in this paper. MHD simulations of accretion discs, such as 
the ones used in our work, generically show significant vari- 
ability due to the stochastic nature of MRI-generated tur- 
bulence and magnetic reconnection. Together, these effects 
can sometimes lead to large flares, similar to those observed 
in the solar corona. Since our goal is to depict "represen- 
tative" spectra, we have to make a choice about what to 
consider "typical" behavior. One option would be to pro- 
duce many simulations using the same initial setup, but with 
different random seeds, and then select only those simula- 
tions that give desirable results. This was the procedure in 
|Moscibrodzka et al.| ( |2009[ ). 

We have instead chosen a different approach, in that we 
perform only one simulation for each set of parameters, and 
then present the median spectrum for each simulation as 
the representative one. By choosing the median value of all 
spectra, we minimize the impact of "extreme" episodes, es- 
pecially a few very brief synchrotron/inverse Compton flares 
that appear to be attributable to numerical limitations of 
the simulations, most particularly the forced axisymmetry. 
For error bars, we give the "1-sigma" variation about the 
median, i.e. the limits within which 68% of the spectra fall. 
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For example, if we have 50 individual spectra (corresponding 
to 50 different simulation times), as is typically the case in 
this work, then for each spectral energy bin, we drop the 8 
highest and 8 lowest data points. In fact dropping just the 4 
highest and lowest data points already gives similar results, 
but we use the 1-sigma values throughout. 

In this paper we only include spectra computed over the 
inner 15rG of the simulation domain. We have confirmed, 
via comparison with spectra computed using the entire sim- 
ulation domain, that most of the emission comes from this 
inner region. The emission from the jets in our simulations is 
completely negligible because of their extremely low density. 
This extreme level of magnetic domination is likely partially 
a byproduct of adhering to ideal MHD, where matter cannot 
effectively load the jets. We hope to explore the jet contri- 
bution to the spectrum more in future works. 



3 SIMULATION SETUP 
3.1 Initial State 

We start each simulation with a torus of gas around a com- 
pact object situated at the origin. The mass of the cen- 
tral compact object is set to the mass of Sgr A* (Mbh = 
4.3 x 10 6 Mq), and the initial density profile inside the torus 
is chosen to produce a desired mass accretion rate, some- 
where in the range 10 -9 and 1O _7 M0 yr _1 (depending on 
the simulation) at the inner grid boundary. 

The free parameters that describe the torus are the 
black hole spin a„ = a/Mmi = cJ/GM^ (where J is the 
angular momentum of the black hole and re = GMbh j(? — 
6.35 x 10 n cm ~ 2.06 x 10" 7 pc ), the inner radius of the 
torus (r- m — 15ro), the radius of the pressure maximum 
of the torus (r cen ter = 25ro), and the power-law exponent 
(q = 1.68) used in defining the specific angular momentum 
distribution, 



9t4> +• 



7/ e , 



8/2-1 



(11) 



We then follow the procedure in Chakrabarti ( 1985 1 to 
solve for the initial internal energy distribution of the torus 
e(r, 6). This sets its initial temperature profile, To = (r — 
l)(/im_r//fcs)e. For the purpose of initialization, we assume 
an isentropic equation of state P — pe(T — 1) = K(F , so 
that now the density is given by p — [e(T — Vf/n] 1 ■ 
We can then use the parameter k to set the density (and 
mass) normalization of the initial torus. 

The torus is seeded with weak poloidal magnetic field 



loops to drive the magnetorotational instability (MRI) ( Bal 
bus fc Hawley|p9T| . The non-zero spatial components are 
B r = —dgA,!, and B = d r A^, where 



Ad 



C(p- pent) 2 sin [AN log(r/S)] for p > p cut , 
for p < p cut , 

(12) 

N is the number of field loop centers, and S — l.lri n . In 
this work we consider configurations with N = 1 and 4, as 
illustrated in Figure [l] in order to investigate the effect of 
changing N. The parameter p cut = 0.5p max .o is used to keep 
the field a suitable distance inside the surface of the torus 
initially, where p ma x,o is the initial density maximum within 
the torus. Using the constant C, the field is normalized such 



that initially /3 mag = P/Pb /3mag,o = 10 throughout the 
torus, where Pb is the magnetic pressure. This normaliza- 
tion is to ensure that the initial magnetic field is weak. This 
way of initializing the field is slightly different than what 



other groups have done, e.g. Beckwith et al. ( 2008 1, who used 
a volume integrated /3 mag , or McKinney &i Blandford ( 2009 1 , 



who used /3 mag = -Rwg/ Pb, a v g , to set the field strength (see 
|McKinney et al.|2012| for a discussion of the different meth- 
ods). We also tested one simulation with /3 mag ,o = 50. The 
higher /3 case ends up with a density (over r < 15rG) that 
is about a factor of 2 higher, a temperature that is 30-50% 
colder, a scale-height that is about 40% thinner, and a mag- 
netic pressure that is about 50% smaller. The higher /3 case 
contributes ~70% more flux in the infrared waveband, and 
~86% more in the X-ray band, compared to the weaker case 
(see D12). 

The fact that f3 mag .o affects our final result is not sur- 
prising since our simulations do not reach a saturated (mag- 



netically arrested or magnetically choked) state l 


Igumen- 


shchev et al. 2003; Igumenshchev|2008; Tchckhovskoy et al. 


2011 


McKinney et al.||2012 


1. Our resulting field strength is 



dictated by how much magnetic field is made available to 
the black hole through our initial conditions. 

In the background region not specified by the torus so- 
lution, we set up a low density non-magnetic gas. Numerical 
floors are placed on density and energy density with the fol- 



lowing forms: Pfloor = 10 p n 



and efl oor = pe = 



10 p ma x,o»" _ • These floors are never applied within the 
disc, nor in most of the background region. They are only 
seldom applied very close to the outer boundary (the few 
last external cells), and more frequently along the vertical 
axis in the funnel region. The most commonly applied floor 
is that on the ratio (p + pe)/Ps- Whenever this quantity 
drops below 0.01 (almost always within the jets), both p 
and pe are rescaled by a factor appropriate to maintain this 
ratio. 



3.2 Grid 

All of the simulations are performed in 2.5 spatial dimen- 
sions (all three spatial components of vector quantities are 
evolved in a single azimuthal slice, although symmetry is 
assumed in the azimuthal direction) using a spherical polar 
coordinate grid. The grid used in the majority of the sim- 
ulations consists of 256 radial zones and 256 zones in polar 
angle. 

In the radial direction, we use a logarithmic coordi- 
nate of the form r\ = 1.0 + ln(r/rBH), where tbh = ra(l + 
\/l — a*) is the radius of the black hole horizon; therefore, 
Ar/r — constant. The inner and outer radial boundaries 
are set at 0.9rBH and 120rc respectively. Note that, be- 
cause we use the Kerr-Schild form of the Kerr metric, we 
are able to place the inner radial boundary some number of 
zones (usually 4) inside the black hole horizon. In principle, 
this choice should keep the inner boundary causally discon- 
nected from the flow, since MHD signals cannot physically 
radiate out from within the event horizon. This situation 
is preferable to having the inner boundary of the grid out- 
side the event horizon, which can lead to artificial behavior 
within the flow. The spatial resolution near the black hole 
horizon is Ar m 0.023rG; near the initial pressure maximum 
of the torus, it is Ar w 0.45r<3. We use outflow boundary 
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Figure 1. Initial magnetic field configuration for the four 
poloidal loop case (top panel) and the single poloidal loop case 
(bottom panel). The colour scales with the value of the magnetic 
field. 



conditions at both the inner and outer boundaries (the MHD 
primitive variables are copied from interior zones into ghost 
zones, except the radial velocity component, which is zeroed 
out if it would lead to inflow onto the grid). 

Since we consider a fairly small radial range, there is 
some concern that our jets (discussed in Section 5.1) might 
be affected by reflections of waves off the outer bound- 
ary of the grid. To test this, we performed one simula- 
tion, B4S9T3M9Ce, that used an extended radial grid with 
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Figure 2. Initial setup of the simulation. The colour/saturation 
scales with the density. The number of cells scales inversely with 
distance from the black hole and angle away from the equatorial 
plane. (The left panel shows the full computational domain while 
the right panel only shows the region of interest closer to the 
central SMBH.) 



r ma i = 1.1 x 10 4 rc. We found that most properties of this 
simulation were very similar to the corresponding results 
done on the smaller (default) grid. The deviations were con- 
sistent with those expected from using a different random 
perturbation seed, as was the case here since the new sim- 
ulation used a different processor distribution which affects 
the random seeding. 

In the angular direction, we include the full range 
^ 9 ^ 7r, with reflecting boundary conditions applied at 
the poles (meaning that the MHD primitive variables are 
copied from interior zones into ghost zones, with the sign 
reversed on the 9 component of all vector quantities, plus 
the staggered magnetic field component zeroed out at the 
pole). We use a concentrated latitude coordinate X2 of the 
form 9 = X2 + |(1 — h) sin(2x2) with h — 0.3, which gives us 
a better resolution near the midplane (r cen t cr A9 — O.lrc). 
Figure [2] shows the initial state of the simulation together 
with the grid resolution. 

Our choice of resolution is sufficient to ensure that the 
fastest growing poloidal field MRI modes are well resolved 
[Amri = 2-irvAz/Q > 10Az ( |Hawley et al.|2011| >] for at least 
the first few orbits of the simulation. We also confirm tha t 
amag = ~B r B 4 ' /P B « 0.3, as expected (|Hawley et al.|2011|. 



We also performed select simulations at one- half and at dou- 
ble our default resolution to directly test the numerical con- 
vergence of our results (simulations B4S91 and B4S9h in 
Table 1). We find very little variation between our default 
resolution (such as in simulation B4S9) and its higher reso- 
lution counterpart (B4S9h), suggesting our results are well 
converged. 

Since our ultimate goal is to compare spectra from our 
simulations with actual data, our true criterion for demon- 
strating reasonable convergence is to compare spectra gen- 
erated from simulations at different resolutions. Figure [3] 
demonstrates that simulation B4S9, with a resolution of 
256 x 256, gives a spectrum that agrees with simulation 
B4S9h, with a resolution of 384 x 384, to within our error 
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Figure 3. Median broadband spectra taken from all timesteps 
in the interval 2.5-3.5 orbits for 3 simulations done at different 
resolutions. The error bars represent the 1 sigma limits about 
the median (see text for more details). The spectra contain two 
components: the synchrotron component responsible for the first 
bump, and the inverse Compton component at higher frequencies. 



bars. On the other hand, simulation B4S91, with a resolution 
of 192 x 128, produces a markedly different spectrum that 
diverges from the other two. Recently, |Shiokawa et al.| ( |2012[ ) 
confirmed that once a certain resolution is achieved, further 
increases in resolution have little effect on the spectra. 



3.3 Time interval 

Because the simulations are performed assuming axial sym- 
metry, the anti-dynamo theorem stating that an axisym- 
metric magnetic field cannot be maintained via dynamo 



action (Cowling 19331, prevents the MRI from being self- 
sustaining. After an initial phase of vigorous growth of the 
MRI channel modes, the turbulence gradually dies out as 
the simulations progress. This 2D approximation is a main 
caveat of the simulations we have run, and it is therefore 
important for us to carefully choose the time interval that 
we wish to analyse. We want it to be after the initial phase 
of MRI growth, but before the mass accretion has begun 
to decay too dramatically (longer simulation times are not 
necessarily beneficial in 2D for this reason). Since we have a 
target mass accretion rate in mind for each simulation, we 
can use that as a guide for selecting our time interval. 

Figure [4] shows the variation of the accretion rate as 
a function of time for three simulations that include cool- 
ing (B4S9T3M9C, B4S9T3M8C, and B4S9T3M7C). We see 
that it takes roughly 1.5 orbits to have the accretion reach 
its peak value, and that the accretion dies out (returns 
to the background rate) after about orbit 5, where we 
are referring to the circular orbital period at r = r ccnt cr, 
i.e. t orh = 2ir(g t4j £ - 94,4,) / (gtd - gt$) = 1.67 x 10 4 s. 
For these three simulations, the target accretion rates were 
10" 9 , 10~ 8 , and 6 x 1O" 8 M yr" 1 , respectively. Figure [I] 
then suggests that an appropriate interval would be be- 
tween 2.5 and 3.5 orbits. We therefore use this as the stan- 
dard time interval for analysis, in the rest of this paper, 
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Figure 4. Mass accretion rate at the event horizon of simulations 
B4S9T3M9C, B4S9T3M8C, and B4S9T3M7C as a function of 
time. 



librium over this time interval, based, for instance, on the 
mass flux being constant as a function of radius over this 
region. Therefore, we can also be confident that our selected 
time interval is late enough not to be affected by transient 
solutions. 



4 RESULTS: THE IMPORTANCE OF 
INCLUDING COOLING LOSSES 

As we have stated, our primary goal in this paper is to as- 
sess the importance of including radiative cooling losses self- 
consistently in numerical simulations of LLAGN like Sgr A*. 
We do this by comparing two types of simulations: those 
that include radiative cooling losses self-consistently in the 
dynamics (indicated by a "C" in the simulation name or 
cooling "ON" in Table 1) and those that neglect them (no 
"C" , cooling "OFF" ). We performed a set of 25 different sim- 
ulations to assess the importance of radiative cooling losses, 
as well as study the influence of parameters such as the ini- 
tial magnetic field configuration, the spin, the ion to electron 
temperature ratio, and the mass accretion rate. 

Table 1 summarizes all of the simulations performed. 
The name given to each simulation is derived from its pa- 
rameter choices in the order they appear in the table. We 
also include the average and root mean square (rms) fluctu- 
ations of the mass accretion rates actually measured in our 
simulations. This value is to be compared with our target 
values of M. In simulations without cooling, this value is ar- 
bitrary since the initial density (or mass) of the disc is a free 
parameter (other than the requirement that the total mass 
of the disc Mdisc be negligible compared to Mbh since we 
ignore the self-gravity of the disc). Its choice does not affect 
the subsequent evolution of the simulation, so can be freely 
rescaled after the fact. This freedom does not extend to the 
simulations that include cooling, since the cooling function, 
A, depends directly on p. Therefore, to test different mass 
accretion rates, it is sufficient in simulations without cooling 
to do a single simulation and then simply rescale it to each 
of the target accretion rates, whereas in the simulations with 



© 2012 RAS, MNRAS 000, [T|fl3l 



General relativistic MHD simulations of accretion around Sgr A *: How important are radiative losses? 



Averaged matter density 




***** 



B4S9M9 
B4S9T3M9C 

B4S9M8 
B4S9T3M8C 

B4S9M7 
B4S9T3M7C 



10 15 
Radius [r G ] 



20 



Figure 5. This figure illustrates the relative importance of radia- 
tive losses on disc density. The plot shows the averaged density 
for simulations B4S9T3M9C, B4S9T3M8C, B4S9T3M7C and 
B4S9 (rescaled to three different accretion rates). A time average 
is taken over the interval i m i n — tmax, and at each radius we take 
an average over the shell. 
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Figure 6. Time-averaged matter density (in g/cm 3 ) of the sys- 
tem for simulations B4S9 and B4S9T1M7C, respectively. The 
simulations are initially identical; differences in the images are 
due entirely to cooling losses, included in the simulation in 
the right panel. The time average is taken over the interval 



cooling, a new simulation must be conducted for each new 
mass accretion rate. This physical scaling is an important 
distinction between our work and that of previous authors 
- each of our simulations with cooling has a real, physical 
mass accretion rate associated with it; it is no longer a pa- 
rameter that can be freely fit to the data. When comparing 
simulations with cooling against those without, we always 
assume the same initial accretion rate. 



4.1 Effects of the cooling losses on the dynamics 

Since our cooling function A depends on p, T e , and b 2 (as 
well as H), one way to assess the importance of cooling in 
the simulations is to observe how these quantities compare 
between simulations that include cooling and those that do 
not. Differences would indicate that the simulations with 
cooling are adjusting dynamically to the radiative losses. 
We present these comparisons in Figures [5] - [8] 

Figure [5] shows the time- and shell-averaged density, p, 
for four different simulations: B4S9T3M9C, B4S9T3M8C, 
B4S9T3M7C, and B4S9. The simulation B4S9 has been 
scaled to three different accretion rates to be consistent with 
the three simulations that include cooling. A clear trend 
is apparent in this figure: as the mass accretion rate in- 
creases, the differences between simulations with cooling and 
those without also increases. At our lowest accretion rate 
(M ~ 1O" 9 M yr" 1 ~ 10" 8 A'/ Ed d), the differences are rel- 
atively small (< 30%), while at our highest accretion rate 
(M « 2tt x 10" 8 M Q yr" 1 ~ 6.6 x 10~ 7 Af E dd), they be- 
come more substantial (~ 50 — 70%). Note that the appar- 
ent "bump" in density at r m lbrc is simply indicating that 
the initial mass reservoir (torus) has not fully redistributed 
itself by this time in the simulation (as 7~i n — 15r<3). The en- 
hanced density associated with the cooling simulations can 
also be seen in Figure [6] which compares simulations B4S9 
and B4S9T3M7C. With the cooling losses included, we end 
up with a thinner, denser disc. 

In Figure [7] we show a similar plot for the density- 



weighted magnetic field magnitude, b 2 . Again the trend is 
that the differences between simulations with and without 
cooling become more pronounced at larger accretion rates 
even-though the trend is less obvious. We have a perfect 
match for cooling vs non cooling at the lowest accretion rate 
for the inner part of the disc, while the curves are distinct for 
the higher accretion rates. Also, the overall normalization of 
the magnetic field increases with mass accretion rate, which 
makes sense since the magnetic field is, roughly speaking, 
normalized by the mass density of the fluid. Further, since 
we are assuming ideal MHD, the magnetic field is trapped 
in the fluid, so as the torus spreads into a disc, the magnetic 
flux will be carried with it. Therefore, it is not surprising 
that the differences in Figure [7] are not as large as in Figure 
[5] for p. 

The density-weighted temperature (representative of 
the disc temperature) is shown in Figure [sj which helps il- 
lustrate one of our main conclusions in this paper. The close 
agreement between the non-cooling simulations and simu- 
lation B4S9T3M9C (target M » 10" 9 M© yr" 1 ) suggests 
that at this level of accretion, radiative losses are not impor- 
tant. However, for M > 1O" 8 M yr" 1 i.e. M > 10" 7 M Bdd , 
the differences in temperatures of the discs dramatically in- 
creases, with the highest M simulation almost an order of 
magnitude colder than the lowest, at small radii. This result 
suggests that only for the very lowest luminosity AGN (Sgr 
A* being the only known candidate example at this time) is 
it reasonable to neglect radiative cooling losses in numerical 
simulations. Figure [9] demonstrates the dramatic difference 
in temperature, up to 2 order of magnitude in the region of 
interest (the inner disc), between our highest accretion rate 
simulation, B4S9T3M7C, and a simulation that neglected 
cooling (scaled to the same accretion rate). Note, especially, 
that only with cooling does any gas with T < 3 x fO 11 K 
reach the black hole event horizon. When cooling is turned 
on, the lost radiation pressure results in a thinner disk, and 
this in turn compresses the gas and magnetic fields, so we 
end up with a thinner but denser disk in the inner part of 
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Table 1. Description of simulation parameters 



Simulation 


B loops (TV) 
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ON 
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3 
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ON 
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B4S5T3M7C 


4 


0.5 


3 


6.3 x 10" 8 


2.03±1.85 x 10~ 7 


ON 
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B4S75T3M7C 


4 


0.75 


3 
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ON 
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3 
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3 
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256 x 256 
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6.3 x 10" 8 
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Figure 7. Piot illustrating the relative importance of radiative 
losses on the resulting magnetic field. The plot shows the average 
of Pb for simulations B4S9T3M9C, B4S9T3M8C, B4S9T3M7C 
and B4S9 (rescaled to three different accretion rates). A time 
average is taken over the interval t m i n — t ma x) and at each radius 
we take a density-weighted average over the shell. 



the disk, while the temperature is much cooler due to the 
radiative losses. 



4.2 Effects of the cooling losses on the resulting 
spectra 

Using numerical simulations together with spectral model- 
ing codes to fit Sgr A* has already been considered by many 
authors (e.g. Goldston et al.|2 005 Moscibro dzka et al.|2009 



Dexter et al.||2010| |Hilburn et al.||2010| |Shcherbakov et al.| 



2010). What is new about our work is that it is the first that 



a numerical study of Sgr A* includes the effects of radiative 
cooling self-consistently. 

In our companion paper (D12), we present the details 
of producing broadband spectra from our simulations, to- 
gether with a parameter study that we use to constrain the 
physical environment of Sgr A*. In this paper we focus on 
how including cooling affects the dynamics within the sim- 
ulations more generally, using sample spectra to illustrate 
the results. As in the previous section, we consider three 
different target accretion rates. 

Figure [T7J] shows the median broadband spectral energy 
distribution for four different simulations: B4S9T3M9C, 
B4S9T3M8C, B4S9T3M7C, and B4S9. The simulation B4S9 
has been scaled to three different accretion rates to be con- 
sistent with the three simulations with cooling. The spec- 
tra represent the synchrotron and inverse Compton emis- 
sion of the system at each frequency. They also include all 
special and general relativistic effects via general relativis- 



tic ray tracing using the codes GEOKERR (Dexter & Agol 
[2009| > and grtrans ( |Dexter|[20l7| ). We can see that for the 
low accretion rates (M ~ 10" a M Q /yr = 10" 8 M E dd), the 
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Figure 8. Plot illustrating the relative importance of radiative 
losses on disc temperature. The plot shows the averaged temper- 
ature for simulations B4S9T3M9C, B4S9T3M8C, B4S9T3M7C 
and B4S9 (rescaled to three different accretion rates). A time av- 
erage is taken over the interval t m i n — i maI , and at each radius we 
take a density-weighted average over the shell. Because the tem- 
perature does not change when the density is rescaled, the three 
non-cooling curves (pink, blue, and orange) are exactly on top of 
each other whether M = 10 -9 , 10~ 8 , or 10~ 7 M Q yr _1 
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Figure 9. Time-averaged temperature (in Kelvin) of the system 
for simulations B4S9 and B4S9T1M7C, respectively. The simu- 
lations are initially identical; differences in the images are due 
entirely to cooling losses, included in the simulation in the right 
panel. The time average is taken over the interval t m i n — tmax- 
Note that the system is very dynamical, with many different time 
scales playing themselves out. Even averaging over an orbital pe- 
riod does not give a perfectly symmetric system as analytic mod- 
els would predict. 



spectra of simulations B4S9T3M9C and B4S9 are nearly 
indistinguishable, while for the high accretion rates (M > 
lO _8 M0/yr = W~ 7 MEdd), the spectra obtained without in- 
cluding the cooling losses self-consistently in the simulation 
are quite different than the ones where they are included. 
Thus, in the regime M > KT 7 M E dd, we conclude that sim- 
ulations which do not take cooling losses into account are 
likely not providing a realistic representation of the physics 
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Figure 10. Broadband spectra for simulations B4S9T3M9C, 
B4S9T3M8C, B4S9T3M7C and B4S9 (rescaled to three differ- 
ent accretion rates) with an inclination angle of 85 degrees. The 
median is taken over the interval t m j n — i max , errors as described 
for Fig. [J 
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Figure 11. Broadband spectra for simulations B4S9T3M9C and 
B4S9 (rescaled to fit the point 230 GHz, 3 Jy) with an inclina- 
tion angle of 85 degrees. The median is taken over the interval 
imin — £max- Observational data for Sgr A* are represented in 
red. fFalcke et al.\B000f/ . 



occurring during the accreting process. In this case, fitting 
for parameters such as spin, temperature ratio, magnetic 
filed configuration, inclination, and accretion rate is likely 
to yield incorrect results. This would seem to apply to all 
known LLAGN except perhaps Sgr A*. 

Of the different simulations considered in this work, the 
spectra of simulation B4S9T3M9C is the most compatible 
with the observational data of Sgr A* (see Figure 11 1. Note 



that we are not trying to fit the lowest energy data, as the 
radio emission is believed to be emitted from a region much 
further out than the radial extent of our simulation domain. 
In contrast, the so-called "sub-mm bump" (see, e.g., |Melia| 



fc Falcke]|2001 1 is expected to originate from quite close to 
the black hole, i.e. the region represented by our simulation. 

To provide concrete values allowing comparison with 
other works, we find that we can fit the Sgr A* data at 230 
GHz, which has a value of 3 Jy, or vL v ~ 4 x 10 35 erg/s, with 
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an accretion rate of 2.60 ± 1.54 x 10 9 M0/yr via the non- 
scalable, cooled simulation B4S9T3M9C (see Table 1). To 
fit this data with the corresponding non-cooling simulation 
B4S9, we must scale the simulation to have an accretion rate 
of 2.95 ±0.87 x 10" 9 M Q /yr at the BH horizon. These values 
are statistically consistent, confirming our earlier point that, 
at the currently best fit accretion rate, it is not necessary to 
consider radiative cooling processes when simulating accre- 
tion onto Sgr A*. For a complete parameter study and the 
best fits to Sgr A* with our models, see D12. 



5 RESULTS: PARAMETER STUDY 

Because we have considered a fairly large parameter space in 
our study, we briefly report how our simulations are affected 
by each of these. Many of the conclusions reached in this 
section confirm results of earlier studies by other authors. 
We include them again here (with references to the earlier 
works) for completeness and to allow comparison. 



Table 2. Jet energies and efficiencies 



5.1 



Influence of the initial magnetic field 
configuration 



As noted in previous studies (e.g. Beckwith et al. 2008 McK- 



mney fc Bla ndford 2009), the initial magnetic field config- 



uration has a major influence on the launching and contin- 
ued powering of jets in numerical simulations. Although our 
study is more focused on the disc than the jets, we did make 
a measure of the instantaneous energy carried in the fluid 
and magnetic components of the jets (defined as material 
with b 2 1 p > 1) at the time of each data writeout from our 
simulations. The fluid and magnetic energy of the jet are 
defined as 



E 



FLjetl, 



and 



mo + P)] drd6d4> (13) 



g (2P B u°u + P B - b°b ) drd9d<t> 



-EEM,jct(i) 

(14) 

Table 2 shows these quantities, averaged over the usual in- 
terval (£ m in — imax), for a subset of our simulations. We also 
report the s ame measure of jet efficiency as Tchekh ovskoy| 
et al.|(|20li|>, namely rj = (M — E)/{M), where 



M{t) 



is the mass accretion rate and 



-gpu d9d<f> 



^/Tld6d4> 



(15) 



(16) 



is the energy flux, both taken at the black hole event hori- 
zon tbh- Angle brackets indicate a time-averaged quantity. 



The negative sign in equation ( 15 I indicates that the mass 
flux is positive when rest mass is transported into the black 
hole. Similarly, E is constructed such that a positive value 
indicates a net flux of energy into the black hole. 

In our models we start with a relatively weak (^ 17. 5G) 
magnetic field contained entirely within the disc, initially in 
the form of poloidal loops. We tested two different configu- 
rations: a single set of poloidal loops centered on the pres- 
sure maximum of the torus and following contours of pres- 
sure/density (an example is simulation B1S9T3M9C) and 



Simulation 



■^FLjet 

(erg) 



^BMJet 

(erg) 



B4S9rT3M9C 


3.43 


X 


10 38 


6.23 


X 


10 38 


0.225 


B4S0T3M9C 


1.40 


X 


10 39 


1.88 


X 


10 39 


0.0625 


B4S5T3M9C 


7.23 


X 


10 39 


9.80 


X 


10 39 


0.120 


B4S7T3M9C 


5.50 


X 


10 39 


1.07 


X 


10 40 


0.161 


B4S9T3M9C 


1.03 


X 


10 40 


2.19 


X 


10 40 


0.464 


B1S9T3M9C 


1.07 


X 


10 41 


6.47 


X 


10 41 


0.559 


B4S98T3M9C 


1.72 


X 


10 40 


7.60 


X 


10 41 


0.525 



four sets of poloidal loops spaced radially, with alternating 
field directions in each successive loop (an example is simu- 
lation B4S9T3M9C). 

The 1-loop configuration yields a significantly more en- 
ergetic outflow (mostly associated with the "jets") - more 
than an order of magnitude stronger than for the more com- 
plex, though perhaps more realistic, configuration consisting 
of 4 initial loops. We have confirmed that, as expected, the 
overall power of both components scales roughly linearly 
with the mass accretion rate, meaning the efficiencies of the 
outflows are not strongly effected by the cooling. 

Aside from the outflow, the initial magnetic field con- 
figuration also has a strong influence on the accretion rate 
as seen by comparing the measured M in Table 1 for simu- 
lations B4S9T3M9C and B1S9T3M9C. The 1-loop simula- 
tion has a substantially higher average M, as well as much 
larger statistical fluctuations. This difference is likely a con- 
sequence of the axisymmetric nature of these simulations. 
In our simulations, accretion is driven by the MRI. For the 
1-loop case, the dominant channel solution is able to disrupt 
nearly the entire disc, driving large amounts of material onto 
the black hole in multiple short, intense episodes. In the 4- 
loop cases, the channel solutions are restricted in their radial 
range and can not disrupt the disc as effectively; accretion 
proceeds more smoothly, and at an overall lower normal- 



ization. [The channel solution ( Hawley fc Balbus||1992 1 is a 
particularly violent form of the poloidal MRI and is char- 
acterized by prominent radially extended features.] These 
effects are illustrated in Figure |12"| 

Figure [13] compares the spectra for the two different 
initial magnetic field configurations, the "1-loop" configura- 
tion (in black) , and the "4-loop" configuration (in red) . It is 
clear that the magnetic field configuration has a strong influ- 
ence on the resulting spectra and is therefore one of the big 
uncertainties affecting all such models of the accretion flow. 
The lack of a priori knowledge of the field configuration 
introduces an effective error in all simulation conclusions 
that is comparable to the difference introduced by including 
cooling. Thus, while we can treat the cooling losses self- 
consistently, it is difficult to overcome the variations that 
different choices in initial magnetic field configurations will 
introduce. 



5.2 Influence of the Spin 

Via its formation and accretion history, Sgr A* is likely a 
rotating SMBH, but as with all black holes it is currently dif- 
ficult to reliably determine its spin. Although the observed 
spectrum should depend on spin, it also depends on the mass 
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Figure 12. Density (in g/cm 3 ) of the system for simulations 
B1S9T3M9C and B4S9T3M9C at the same time step (t= 3 or- 
bits). The two simulations started with the same initial torus 
configuration as shown in figure [2] The only initial difference 
between the two models is that the left panel started with a sin- 
gle magnetic loop configuration (as shown in the bottom panel 
of /igureQp, while the right panel started with a 4 magnetic loop 
configuration (as shown in the top panel of figure^. We can see 
that in the case of the single initial magnetic loop, the channel so- 
lution extends through the entire disc, resulting in a higher mass 
accretion rate at the black hole horizon and a more extended disk. 
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Figure 13. Broadband spectra for simulations B4S9T3M9C, 
B4S9T3M9, B1S9T3M9C, and B1S9T3M9. The median is taken 
over the interval t m i n — t m ax- 



accretion rate, which is somewhat constrained, as well as the 
geometry of the accretion flow. Nevertheless, spectral fitting 
of Sgr A* has already placed some meaningful constraints 
on the value of a* ( Broderick et al.|2011 Shcherbakov et al 
|2010[ ), although there can be additional complicating fac- 
tors, such as disk tilt ( [Dexter fc Fragile|[20TI| |2012| . Other 
methods of inferring the ISCO, and then the spin, such as 
observing relativistically-broadened line profiles or modeling 
the thermal emission from the disc, require the presence of a 
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Figure 14. Jet fluid energy as a function of spin. The six 
data points represent the averaged jet energies for B4S9rT3M9C, 
B4S0T3M9C, B4S5T3M9C, B4S7T3M9C, B4S9T3M9C, and 
B4S98T3M9C. 
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Figure 15. Jet magnetic energy as a function of spin. The six 
data points represent the averaged jet energies for B4S9rT3M9C, 
B4S0T3M9C, B4S5T3M9C, B4S7T3M9C, B4S9T3M9C, and 
B4S98T3M9C. 



geometrically thin, optically thick disc, which is not observed 
at the low accretion rates of Sgr A*. The ultimate determi- 
nation of spin will likely require future observations from 
mm-VLBI (poeleman et al.|[2008] |Fish et al.|[20TT] |Broder-| 



ick et al.|2011[ ). 

We investigated the spin dependence of our results by 
performing simulations at six different values of a*: -0.9, 0, 

0. 5, 0.75, 0.9, and 0.98, where the minus sign indicates ret- 
rograde rotation (disc and black hole rotating in opposite 
senses). We will discuss the effects of spin on our spectral 
fits in companion paper D12. Concerning the disc evolution, 
there does not seem to be any clear trend in M, see Table 

1. If anything there appears to be a peak in the distribu- 
tions, with the highest mass accretion rates occurring at 
intermediate spins, 0.3 < a, < 0.7. The outflow power on 
the other hand, shown in Table 2, reveals a clear trend with 



spin, as expected based on theoretical grounds (e.g. Bland- 
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ford & Znajek 1977 Tchekhovskoy et al. 


2010 


1 and Irom 


previous numerical work (De Villiers et al. 


2005 


Hawley & 


Krolik 2006, McKinney 2006; Tchekhovskoy & McKinncy 



2012). This trend is seen in both the fluid and magnetic 
components of the outflow and the correlations are shown 
in Figure [M] 

The effect of spin on jet power is a question that has 
been debated recently in the literature, specifically using 
claimed spin observations from black hole X-ray binaries. 
So far the observational results are not definitive (see, e.g., 
the opposing conclusions found in |Fender et al.|[2010| and 
Narayan & McClintock 20121. Our results, while focusing 



on SMBHs, would be predicted to scale to smaller mass 
black hole systems via the Fundamental Plane of BH accre- 
tion relation ( |Merloni et al.|2003| [Falcke et al.|2004[ ). Semi- 
analytical models for the compact jets observed in weakly 
accreting systems predict a dependan ce of radio luminos 
ity on total jet power of ~ Q 1 ' 4 (e.g 



Falcke & Biermann 



1995 ). Thus if we consider the total power as the sum of the 
fluid and magnetic components in Figure (15} we would also 
predict a correlation between observed radiative power and 
spin. 



5.3 Disk Thickness 

In all our simulations we started with the same disk thick- 
ness (Ar/r =constante ~ 1) in order to start with a thick 
disk, appropriate for the Galactic center. Then the shape 
of the disk is modified during the simulation and is espe- 
cially sensitive to the cooling as shown in figure [6] We also 
investigated thinner disk simulation, the thinner disk em- 
mit a little bit less that the thick disk but they could be 
compatible within one sigma. 



6 CONCLUSIONS 

For the first time we have been able to assess the importance 
of the radiative losses in numerical simulations of LLAGN, 
specifically Sgr A* [see also Moscibrodzka et al. (2011 1], We 



show that radiative losses can affect the dynamics of the 
system and that their importance increases with accretion 
rate. We set a rough limit of M > 10 _7 AfEdd, above which 
radiative cooling losses should be included self-consistently 
in numerical simulations. Otherwise, many important de- 
rived dynamical quantities, such as density, magnetic field 
magnitude, and temperature, may be off by an order of mag- 
nitude or more, especially when the accretion rate reaches 
M ~ 10" 6 M E dd, which correspond to 10" 7 M Q yr" 1 for Sgr 
A*. Since several recent works suggest that accretion physics 
is similar across the mass scale, this accretion rate in Edding- 
ton units should likely be an important limit for all black 
holes. Thus, we predict that the inclusion of self-consistent 
radiative cooling above ~ 10 _6 MEdd should be important 
for LLAGN in general. 

Overall, this study allows us to have a more consis- 
tent model of accretion, even for the well-studied and under- 
luminous source, Sgr A*. Not only do we have a more real- 
istic model by including the physics of radiative losses, but 
we are also able to show the influence of the accretion rate 
on the resulting spectra. 



The spin of the central black hole and the initial mag- 
netic field configuration of the torus also have important 
consequences on the dynamics of the system and the result- 
ing spectra. By including the cooling losses in our study, we 
can discuss the influence of these free parameters with more 
accuracy. For instance, initial magnetic field configurations 
consisting of a single set of poloidal loops result in signifi- 
cantly more powerful outflows than the four-loop cases, and 
we find that the jet power increases with the spin of the 
central black hole. 

As mentioned in Section |3.3| there is concern as to 
whether or not our simulations have run long enough to 
reach meaningful equilibrium states. This concern is espe- 
cially pertinent for our case using 2.5D simulations, as these 
can never truly reach a steady state. The reason is that, after 
a period of initial vigorous growth, the MRI in 2.5D simu- 
lations steadily decays because the dynamo action that nor- 
mally sustains it requires access to non-axisymmetric modes 
that are obviously inaccessible. We, therefore, chose a time 
interval for analysis when the mass accretion rate closely 
approximated our target value for most simulations. Clearly 
simulations in 3D with a longer duration (to ensure a proper 
equilibrium is reached) would be beneficial, for comparison. 

One other concern with only performing 2.5D, axisym- 
metric simulations is that it precludes the possibility of ex- 
ploring the effects of misalignment between the angular mo- 
mentum of the gas and the black hole. In reality, most super- 
massive black holes are unlikely to be accreting from matter 
sharing the same orbital plane as the black hole spin. |Dexter| 
& Fragile (20121 recently showed that accounting for such 



a "tilted" accretion flow can dramatically alter the best-fit 
characteristics of Sgr A* and produce important new fea- 
tures in the spectra. 

As a final note, one can justify neglecting full radia- 
tive transfer in simulations of Sgr A* and likely most other 
weakly accreting LLAGN because the inner regions are gen- 
erally thought to be optically thin. However, to treat the 
outer regions of the accretion flow, or higher luminosity 
sources, a more thorough treatment of radiative transfer will 
need to be implemented into the simulations. Similarly, to 
approach important questions such as the mass loading and 
particle acceleration in the jets, likely resistive (non-ideal) 
MHD will need to be considered. We see this work as an 
important step towards these next technological "horizons" . 
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